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We study a reweighting method aiming at numerical studies of QCD at finite density, in which 
the conventional Monte-Carlo method cannot be applied directly. One of the most important 
problems in the reweighting method is the overlap problem. To solve it, we propose to perform 
simulations at several simulation points and combine their results in the data analyses. In this 
report, we introduce this multipoint reweighting method and test if the method works well by 
measuring histograms of physical quantities. Using this method, we calculate the meson masses 
as continuous functions of the gauge coupling j3 and the hopping parameters K in QCD at zero 
density. We then determine lines of constant physics in the ()3, K") space and evaluate the deriva¬ 
tives of the lattice spacing with respect to j3 and K along the lines of constant physics (inverse of 
the beta functions), which are needed in a calculation of the equation of state. 
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1. Introduction 

In a study of QCD phase diagram at finite temperature (T) and density, the complex quark 
determinant causes a serious problem in numerical simulations. The reweighting method is com¬ 
monly used to avoid this problem in the low density region. However, when we increase the chemi¬ 
cal potential (/i), the sign problem and the overlap problem becomes severe. In this report, we focus 
on the overlap problem. The overlap problem is expected to be milder if one changes a couple of 
parameters at the same time. For example, WHOT-QCD collaboration investigated the phase struc¬ 
ture of A/^f-flavor QCD in the heavy-quark region, and found that the system at large quark masses 
is controlled by only two combinations of parameters, j3 Kj and Y!fLi cosh(jU//r) 

on an xNt lattice, where jS = 6/g^, and Kf and jXf we. the hopping parameter and chemical 
potential of /^-flavor, respectively [p. This means that, when one changes the parameters with 
keeping these combinations constant, the system does not change, thus the overlap problem does 
not arise. We expect such combinations of parameters exist also in the light-quark region. 

In that study in the heavy-quark region |[I]], the multipoint reweighting method [§] for j3 played 
an important role: Combining data obtained at different j3, we could calculate the effective potential 
in a wide range of the observable values, that was mandatory in a reliable evaluation of the transition 
point. In the present paper, we extend the multipoint reweighting method in the multi-parameter 
space of j8 and Kf to overcome the overlap problem in the light-quark region. 

In the next section, the multipoint reweighting method is introduced, and, in Sec.^, we examine 
the overlap problem performing numerical simulations in two-flavor QCD. We then calculate the 
meson masses, the lines of constant physics, and the derivatives of the lattice spacing with respect 
to j8 and k along the lines of constant physics in Sec. The section ^ is for the conclusions. 

2. Multipoint reweighting method 

We extend the multipoint reweighting method to the multi-parameter space such that both 
gauge and quark parameters are changed simultaneously. Let us first define a histogram for a set of 
physical quantities X = (Xi ,^ 2 , ■■■) as 

r Nf 

detM(fC/,/iy). (2.1) 

4 i f=i 

where Sg is the gauge action, M is the quark matrix and X = {Xi,X 2 , ■ ■■) is the operators for X. 
We denote S = Sg — ^ylndetM. The coupling parameters of the theory are jS, fc = (fCi, k: 2 , • • •, Kjvj), 
and jU = (jUi,/X 25 • • • iMWf)- Then, the partition function is given by Z(j3, fc,/l) = / w{X-,l5,K,jl)dX 
with dX = dXi, and the probability distribution function of X is given by ^w(X; jS, ic, /I). The 
expectation value of an operator ^[X] which is written in terms of X is calculated by 

mx])ip.K,n) = zili^K,^) I ( 2 . 2 ) 

For simplicity, we denote the set of coupling parameters (j3, fc, /I) as jS in the following of this 
section. 
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We perform a simulation at j3o and calculate the histogram at j3. We fix three quantities X, 
S{[5) = S and S{po) = So to construct the histogram, where S is the value of the action with the 
coupling parameters j8 on the configuration generated at jSo- From Eq.(^^, we find 

w(X,S,So;i3) =e-(^-^«)w(X,S,5o;iSo). (2.3) 


The histogram ofX is given by w(X;j8) = / w{X,S,So‘,l5)dSdSo. 

When X has large correlafion wifh fhe difference of fhe actions, S — Sq, the peak position of 
the distribution may change appreciably, causing the overlap problem. (See Sec. ^.) To overcome 
the overlap problem and to obtain w and Veff that are reliable in a wide range of X, we extend 
the reweighing formulas to combine data obtained at different simulation points [|^, ^ for the case 
including the fermion action. 

We combine a set of Ngp simulations performed at j3, with the number of configurations Ni 
where / = 1, • • • ,Nsp. Using Eq. (^), the probability distribution function at pi is related to that at 
P as 


Z-\pi)w{X,S,S-pi)=Z-\pi)e-^^‘-^^w{X,S,lp). 


(2.4) 


Summing up these probability distribution functions with the weight Ni, 




No 


'£NiZ-^iPi)w{X,S,lpi)=e^'£NiZ-\pi)e-^-w{X,S,lp), 

we obtain 

wiX,S,lp) = G{S,S-,p,p)Y,NiZ-\pi)wiX,S,lpi[ 
where ^ = {Pi,---,PnJ and 


i=l 




1=1 


(2.5) 


( 2 . 6 ) 


G{S,S-p,p) = — -. 


(2.7) 


Note that the left-hand side of Eq. (2.5) gives a naive histogram using all the configurations disre¬ 
garding the difference in the simulation parameter (j3, if). The histogram w{X,S,S',p) at (j3, if) is 
given by multiplying G{S,S-,p,p) to the naive histogram. 

The partition function is given by 






Z{P) = '£Ni GiS,S-p,p)Z-\pi)w{X,S,S-pi)dXdSdS='£Ni{G{S,S-p,p)) . (2.8) 


1=1 


1=1 


(A) 


The right-hand side is just the naive sum of G{S,S',p,p) obtained on all the configurations. The 
partition function at Pi can be determined, up to an overall factor, by the consistency relations. 


No 




-Si 


Z{Pi) = Y^NUG{S,S-Pi,P)) = ^ 

ti ^ Gh) tx yLdNje-^>Z-\pj)^ 


(2.9) 


(ft) 
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for / = 1, • • • jA^sp- Denoting ft = — lnZ(j3,), these equations can be rewritten by 

1 = Y!k^iNk I — Sj — fi + fj]\ \ . Starting from appropriate initial values of 

\ ^ ^ / (ft ) 

fi, we solve these equations numerically by an iterative method. Note that, in this calculation, one 
of fi must be fixed to remove the ambiguity corresponding to the undetermined overall factor. 

The expectation value of an operator X at f can be evaluated as 

N- 

{X)ip) = ^Jxw{X,S,S-li)dSdS = ^^f^Ni(xG{S,S-l5,f))^^^^, (2.10) 

and the histogram of X is obtained by 

w(X;j8) = Y^Ni{5{X-X)G{S,S-15,15)) (2.11) 

Note that, again, ) the right-hand side is just the naive sum of XG over all the 

configurations disregarding the difference in the simulation point j3,. 

3. Overlap problem and histograms 

To study if the multipoint reweighting method is useful in avoiding the overlap problem, we 
perform simulations of QCD with degenerate 2 flavors of clover-improved Wilson quarks coupled 
with RG-improved Iwasaki glue, at /r = 0. The improvement parameters of the action are the same 
as those adopted in Ref. 0. The simulations are carried out on an 8^ lattice at 9 simulation points 
(3 j8’s X 3 fc’s) for the test study in this section, and on a 16^ lattice at 30 points (6 jS’s x 5 tc’s) in 
Sec. ^ The number of configurations for the measurement is 200 at each simulation point. 

The action S required in the reweighting method is composed by Wilson loops and IndetM. 
The calculation of the fermion part requires large computational cost. In this study, we evaluate 
detM by measuring the first and second derivatives of IndetM at several tc^’s, where i = 1,2, • • •, 
on each configuration, and interpolate the determinant between fc, and assuming a quadratic 
function, IndetM(K') = lndetM(K',)-|-C i(k'— fc,)-|-C 3 (ff—fc,)^-|-C 3 (k'— fc,)^-|-C 4 (k'— fc,)'*. The 
coefficients Ca are determined such that the first and second derivatives are consistent with the 
measured values at fc, and Ki^i. The clover term with the coefficient cgw is also evaluated in terms 
of its derivatives. Here, because cgw depends on jS in our choice, it affects the j3-dependence of 
observables. However, we find that the effects from cgw are very small in the range of j3 rele¬ 
vant in our study. We thus approximate the j3-dependence of the action by a linear function, i.e. 
lndetM(j8, k) = IndetM^fo, k) + — po){dcsw/dl5){d IndetM/dcsw)ipQ^'i^) [ll- Note that the 

overall constant of IndetM is not needed. 

In the left panel of Fig. [I[ we show the results for the improved plaquette P = -|- 

2ciVT'^^ of the Iwasaki action at (5 = 1.825, where ci = —0.331, cq = 1 — 8ci, and is the 
(/ X j) Wilson loop. Black dots represent the expectation values of P at the thr'ee simulation points 
without reweighting. Blue, green and purple curves are the results of the naive reweighting method 
using the data at fc = 0.1400, 0.1425, and 0.1440, respectively. Each result of the naive reweighting 
method is reliable around the corresponding simulation point, but fails reproducing far away sim¬ 
ulation results. The reason can be easily understood by consulting the histogram of P: Red curve 
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Figure 1: Left; The expectation value of P = + 2c\W^^^ as a function of tc at j3 = 1.825. Middle: 

The histogram of P at various if’s obtained by the naive reweighting method using the configurations at 
K = 0. 140. Right: The histogram of P by the multipoint reweighting method. 


(P,k)=(1 .8325,0.14000) - {P,k)=(1 .8250,0.14000) 

(P,k)=(1 .8400,0.14000) {P,k)=(1 .8250,0.14160) 

(P,k)=(1 .8475,0.14000) - {P,k)=(1 .8250,0.14240) 




Figure 2: The j3-dependence (left) and ff-dependence (right) of the histogram for and 

N^^^[dS/dK]svB=Nf^i[dS/dK-i288NfKyco)idS/dl5)], where iVsite = 


in the middle panel of Fig. [I] is the original histogram at (j8, tc) = (1.825,0.140), and green, blue, 
magenta and light blue curves are the histograms at fc = 0.1412, 0.1424, 0.1436, and 0.144, respec¬ 
tively, estimated by the naive reweighting method Eq. ( |2.3[ ) using the data at K = 0.140. Because 
the histograms at K other than the simulation point is calculated as the product of the reweighting 
factor and the original histogram, the histograms are not reliable out of the range of the original 
distribution. In fact, the value of P distributes between 1.64 and 1.69. Even when K is changed by 
the reweighting method, the upper and lower bounds of the distribution does not change in Eig. ^ 
(middle). Since the expectation value is approximately the peak position of the histogram, the 
expectation value also cannot go out of the range of the distribution, as shown in Eig. |I| (left). 

To enlarge the range of the distribution, we combine the simulation data obtained at fC = 
0.1400, 0.1425 and 0.1440 using the multipoint reweighting method explained in the previous 
section. The red curve in the left panel of Eig. [I] is the result of the multipoint reweighting method. 
We find that the red curve smoothly connects all the direct simulation results with small error bars. 
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Figure3: Left: The pseudoscalar and vector meson mass ratio mps/my. Right: The lines of constant physics 
for each nips /in the (j3, k) plane. 


In the right panel of Fig. |I|, histograms from the multipoint reweighting method are plotted for 
k: = 0.1412, 0.1424, 0.1436, and 0.144. 

The method is applicable to other observables. The expectation values of dS/d[5 and 
[dS/dK]svB^dS/dK-{28m{Kyco)idS/dl5) 

are needed in the calculation of the equation of state in the integral method. In the right and 
left panels of Fig. ^ we show the j8 - and tc-dependence of the 2-dimensional histogram of these 
quantities using the multipoint reweighting method. The histogram moves as j8 and fc are varied. 
We can compute the expectation values of dSjdp and [^S/^tcJsuB as continuous functions of jS 
and K without the overlap problem. Once we obtain physical quantities as continuous functions of 
coupling parameters, we can calculate the lines of constant physics in the coupling parameter space 
as well as the beta functions. 

4. Lines of constant physics and beta functions 

In this study, we define the lines of constant physics by fixing the dimension-less ratio of 
pseudoscalar and vector meson masses nips/niv in the (j8, tc) space. On the lattice, we measure 
dimension-less observables mpsa and niya, where a is the lattice spacing. Because these lattice 
observables vary as we change (3 or K, a is varied when we move along a line of constant physics. 
The beta functions dp jda and dxjda are defined through the variation a along a line of constant 
physics. The beta functions are needed in the calculation of the equation of state. 

The multipoint reweighting method is useful for the calculation of the beta functions because 
we can calculate observables as continuous functions of j3 and fC. Combining the data at 30 simu¬ 
lation points (6 j3’s X 5 fc’s) on the 16^ lattice by the multipoint reweighting method, we obtain the 
mass ratio mps/my plotted in the left panel of Fig. From these data, determine the lines of con¬ 
stant physics for mps/niy = 0.70, 0.72, 0.74 and 0.76, as shown in Fig. ^ (right). We then calculate 
the derivatives, {mya)dp/d{m\a) and {mya)dK/^{m\a)^kmg each line of constant physics. The 
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Figure 4: The beta functions: {m\;a)d^ / d{mya) (left) and {mya)dK/d(m\/a) (right). 


results are shown in Fig. Combining these beta functions with the measurement of {dS/dx) and 
{dS/dj5) on finite-temperature lattice, we can calculate the equation of state. 


5. Conclusions and outlook 

We discussed the multipoint reweighting method in a multi-dimensional parameter space to 
avoid the overlap problem. Using the method, we can reliably calculate histograms of physical 
quantities as well as the expectation values of the physical quantities, as continuous functions 
of coupling parameters. These enable us to compute the lines of constant physics and the beta 
functions, which are needed in a calculation of the equation of state. 

Our final objective is a study of finite density QCD. Using the multipoint reweighing method, 
we may absorb the main effect of the chemical potential by a change of j3 and K. If so, we may 
investigate finite density QCD avoiding the sign problem. 
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